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Abstract 

We describe methods for computing the Kummer function U{a,b,z) 
for small values of z, with special attention to small values of b. For these 
values of b the connection formula that represents U{a,b,z) as a linear 
combination of two lU-functions needs a limiting procedure. We use the 
power series of the iFi-functions and consider the terms for which this 
limiting procedure is needed. We give recursion relations for higher terms 
in the expansion, and we consider the derivative U'(a,b,z) as well. We 
also discuss the performance for small \z\ of an asymptotic approximation 
of the Kummer function in terms of modified Bessel functions. 
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1 Introduction 

As is well known, confluent hypergeometric functions appear in a vast number 
of applications in physics and engineering. In spite of their importance, few 
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algorithms are available for the computation of any of the standard solutions 
of the Kummer equation in the case of real or complex parameters; see for 
example [71IH] for the computation of M(a,b,x). In this paper we describe a 
method for computing the Kummer function U{a,b^z) for small values of \z\, 
|a| and |&| by using a simple expansion of a function in terms of reciprocal 
gamma functions. To compute the Kummer function t/(a, b, z) for small values 
of l^l, the connection formulas for this function in terms of the -functions 
is the most obvious starting point. When also small values of the ^-parameter 
are allowed, special care is needed because certain terms in these formulas will 
become singular, although the limit for 6—^0 will exist. In case b — 0 the power 
series expansion of U{a, b, z) will contain a logarithmic term. This phenomenon 
occurs as well for positive integer values of b, but when the case for 6 = 0 is 
under control, we can use stable recursions m to handle b = 1,2,3,.... It 
should also be noted that, for numerical evaluations with fixed precision, the 
power series expansions can only be used for rather small positive values of z, 
say, for 0 < z < 1. This is due to the behaviour of the ifi-functions, which 
become exponentially large as z —> -l-oo, whereas the [/-function is only of 
algebraic growth in that case. 

With the coefficients provided by the expansion in terms of reciprocal gamma 
functions, we can safely compute the function U{a,b,z) and its z-derivative in 
double precision for 

0<k|<l, |a|<i, |6|<i. (1.1) 

This computational method complements others described in m. where the 
computation of the confluent hypergeometric function U (a, &, z) is considered in 
other (complementary) parameter domains. 

In the paper we also consider an asymptotic approximation of the Kummer 
function U{a,b,z) in terms of modified KT-Bessel functions and discuss its nu¬ 
merical performance for small |z|. An analysis of the stability of the recurrence 
relations satisfied by the coefficients of the expansion is provided as well as few 
tests illustrating its behaviour. 


2 Power series expansions for small \z\ 

We use the representation (see O Eqn. 13.2.E41]) 

zi-'>r(6-l) ^ fa-b+1 


U{a,b,z)= 


T{a-b+l)^ \b 


I z ] + 


r{a) 


-iFi 


2-b 


;z (2.1) 


and the form 


^ f iFiia;b]z) _^_f,iFi{a - b + 1-2 - b; z)\ 

"<“''’'"> = S^(r(i,)r(a-i, + i)“" -r(a)r(2-t)-)' 


These connection formulas are not defined for integer values of 6, although the 
limit exists for, say, 6—^0. 
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We assume z ^ 0 and b ^ Z. We expand the iFi-functions and obtain 


U{a,b, z) = 


r(i - b) 
r(a —5+1) 


nbz 

sin(7r5)r(a)r(a 


T-:^ Wm —. 

0 + 1 ) m! 

' m—0 


(2.3) 


where 

'^m — 5 

'^m 

— (-^771 

Ajn = T{m + 1)T{2 - b + m)T{a + m + 1), (2.4) 

Bm = z~^r(a — 5 + 1 + TO)r(5 + 1 + m)T{m + 2), 

Vm = r(m + 2)r(5 + TO + l)r(2 — 5 + to). 

The Um are well defined in the limit 5 —>■ 0, but for numerical computations 
we need a stable representation of Wm- We derive such a representation that we 
can use for values of a and b in the interval For other a-values we need 

a few recursion steps when this parameter is not too large. For recursion with 
respect to 5, we refer to ^2.31 

We concentrate on wq] other values Wm follow from stable recursions when 
5 is small. Namely, it is easily verified that 


'^m+1 




dm — is^m 


(2.5) 


where 


dm = Cm — {m^ + (a + 2 )to + a + 1)5, 


5m = Cm + (m + 2)(a - 5)5, 

Cm = (m + 1 )(to + 2)(to + 0+1), 


( 2 . 6 ) 


dm = —{m? + 2m{a + 1) + 3a + 1) + (to + 2)5. 

Theorem 4.17 of [5] can be used to check that {um} is a dominant solution of 
(12.5p . which means that when we have a stable representation of Um^ we can 
compute Um+i in a stable way as 5 —>■ 0. The recursion for Vm is straightforward, 
with starting value 


i;o = r(5+l)r(2-5) 


(1 — 5) 7r5 
sin(7r5) 


(2.7) 


2.1 Computing the derivative 

For the derivative of t/(a, 5, z) we obtain from (12.31) 


t/'(a,5,z) 


_ ^ ^ d ^ 

sin(7r5)r(a)r(a — 5 + 1) ’"to! 

^ ^ ^ ^ ^ ' m—0 


( 2 . 8 ) 
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where 


(2.9) 


= (m + l)Wm + zw'^. 

For w!^ we need, using ( 1 ^ . 

= -Bm, ( 2 . 10 ) 

z 

for which we do not need further analysis. 


2.2 Computing the first term Wq 

Following the ideas described in m, we introduce the function 

(r(a + i + 6) " r(a + i)) ■ 

The reciprocal gamma function is an entire function, and we have the expansion 

- OO 

fpi = 1^1 < 

k=l 

where ci = 1, C 2 = 7 , and the rest of the coefficients satisfy (see [U §5.7]) 

(A:-l)cfc = 7Cfe_i-C(2)cfe_2 + C(3)cfc-3-■ ■ ■ + (-l)^C(^-l)ci, A: > 3, (2.13) 

where 7 is Euler’s constant and (C(z) denotes the ((-function. See Table [TJ for 
31D values see m- 

For G(a, h) defined in (|2.11|) we have the representation 


LXJ 

G(a, 6 ) = ^Cfc4, dk = ^{{a + bf-^-o^-^), (2.14) 

fc =2 


which we use for a and b in the interval [— 5 ,^], or for complex values with 
modulus not exceeding i. We have d 2 = = 2a + b, and the recurrence 

relation 

dfc+2 = (2a -I-6)(ife+i - a(a -I-5)dfc, fc > 1. (2-15) 

For a few recursion steps we can use the relations 


(a -I- l)(a -I- 6 + l)G(a -I- 1, 6 ) — (a -I- l)G(a, b) — — ———r, 

r(a + 1 ) 

(a -I- 2)(a + b + 2)G(a + 2,b) = (2a -1-5-1- 3)G(a -1-1,5) — G(a, 5). 


In this case. Theorems 4.6 and 4.7 of [S] are inconclusive with respect to the 
existence of minimal solutions for these recurrences for G(a, 5). 

For Wo we have 


bwo 


r(a-fl) __j,r(a-5-hl) 
r(5+l) ^ r(2-5) 


(2.17) 
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Table 1: Coefficients Ck of the expansion in (12.121) . 


Cl 

C2 

C3 

C4 

cs 

ce 

c? 

cs 

cg 

cio 

Cll 

Cl2 

Cl3 

Ci4 


1 . 00000000000000000000 ; 

0.57721566490153286061; 

-0.65587807152025388108; 

-0.4200263503409523553e-l; 

0.16653861138229148950; 

-0.4219773455554433675e-l; 

-0.962197152787697356e-2; 

0.721894324666309954e-2; 

-0.116516759185906511e-2; 

-0.21524167411495097e-3; 

0.12805028238811619e-3; 

-0.2013485478078824e-4; 

-0.125049348214267e-5; 

0.113302723198170e-5; 


CIS = -0.20563384169776e-6; 
ci6 = 0.611609510448e-8; 
ci 7 = 0.500200764447e-8; 
ci8 = -0.118127457049e-8; 
ci 9 = 0.10434267117e-9; 

C 20 = 0.778226344e-ll; 

C 21 = -0.369680562e-ll; 

C 22 = 0.51003703e-12; 

C 23 = -0.20 5 8 3 26e-13; 

C 24 = -0.53481e-14; 

C 25 = 0.12268e-14; 

C 26 = -0.11813e-15; 

C 27 = 0.119e-17; 

C 28 = 0.141e-17; 


We write 


^-0 = - 1 ) 


r(6 + i) r( 2 - 6 ) 


r (2 - b) 


(2.18) 


and we express some of the gamma functions in terms of the function G{a,b). 
We have 

p/ _ r(a -I- 1) _ 

^ “ l-5r(a-bl)G(a,-6)’ 


and obtain 

Wo = 


r(l + 6) = 
r (2 -b) = 


r(a-f 1 ) 


1 


l-b6G(0,6)’ 

1 - b 

l-bG{0,-b)' 


(6-l)(l-5r(a-f l)G(a,-6)) 

(r(a + 1)(1 - 5)(1 + 6G(0, &))G(a, -b) 


l + {b-l)G{0,b)-G(0,-b) + 


.-b 


-{l-bG{0,-b 


When we use an expansion for 


(2.19) 


( 2 . 20 ) 


,-f) 


“ 1 1 e 

- = mz— 


— 6 In 2 


- 1 


5 In; 


( 2 . 21 ) 


for small values of |61nz|, we have obtained a representation for wq in which we 
can allow small values of b. To start the recursion for Um given in (12.51) we use 
uq = Vo Wo with Vo given in (12.7p . 
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In m the function 



|^ r(z + e) 

I r(z) 



is considered for small values of e. In our notation we have 


re(z) = —r(z + e)G{z — 1, e). 


( 2 . 22 ) 


(2.23) 


Skorokhodov introduced this function for the computation of the Gauss hyper¬ 
geometric function 2 ^"i(a, b] c; z), where for certain values of the parameters a, 
b and c numerical problems arise. Because Skorokhodov needs the evaluation 
of rg( 2 ;) for general complex values of z, more details and expansions are con¬ 
sidered in m- For other papers on the computation of the 2 F'i-function (with 
attention to the computation of functions like re(z) and G(a, &)), we refer to [3] 
and [4] 


2.3 Positive integer values of b 

When we have control of the computations for small values of b, we can ob¬ 
tain the functions U{a^b + k,z) without further extra analysis, although the 
connection formulas in (EB and (12.211 cannot be used for & G Z. 

We have the relations 

U{a,b 1, z) = U{a, b, z) — U'{a, b, z), 

(2.24) 

zU'{a, b + l,z) = bU'{a, b, z) — aU{a, b, z), 

and we conclude that, when b is small, to obtain the values on the left sides no 
special limits are needed, once we have U{a, b, z) and its derivative. In addition, 
when a, b and z are positive, we have U'{a,b,z) < 0, and the recursion is 
perfectly stable. For more computational aspects of the recursion with respect 
to b, we refer to |I3) . 


2.4 Computing G{a, b) by quadrature 


Because I/r(z) is an entire function, we can write the function G{a,b) defined 
in (I2.II|) in the form 




dz 


c r(z -I- I) (z — a)(z — a — b)' 


(2.25) 


where C is a closed contour that encircles the points z = a and z = a -I- 5 
clockwise, say a circle around the origin with radius r > max(|a|, |a-|-5|). Using 
such a circle, we can write 


G{a,b) = ^j 


dO 


-TT r(^) (z-a)(z-a-5)’ 


z = re 


(2.26) 


As explained in [5l §5.2.3] and m, the trapezoidal rule gives excellent results 
in a quadrature algorithm. For an efficient code, we can precompute the values 
of l/r(z) in a set of equidistant nodes on the circle with radius r. 
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Table 2: Relative errors in the relation in (12.271) for a = 0.2 and 6 = 10 
k = 1(1)5, and two complex values of z. 


k 

z = -0.5-0.R 

z = 1 + i 

1 

0.15e-15 

0.25e-15 

2 

O.lOe-14 

0.76e-15 

3 

O.lle-14 

0.26e-14 

4 

0.33e-15 

0.23e-14 

5 

O.lOe-14 

O.lle-14 


2.5 A numerical test 

For a numerical test we have verified the relation 

U{a— l^byz) = {a — b + z)U{a,b, z) — zU'{a, &, z) (2.27) 

for a = 0.2 and a few values of z and b = 10“^^, k = 1(1)5. In Table [2] 
we give the relative error in this relation. We have used a Maple code with 
Digits=16 and the coefficients in Table [I] we stop the summation of the series 
in (12.31) and (12.81) when the absolute values of the terms are less than 10“^®. For 
z = —0.5 — O.li about 15 terms are needed, for z = 1 + * about 20. Because of 
the strong relation between the terms in the series, the evaluation of the series 
is done together. 


3 Asymptotic representation in terms of modi¬ 
fied Bessel functions 

Our starting point is the representations for a large given in [m §10.3] 

i(i-h) 


C/(a,&,z) = 2(0' 


r(a) 


RTh-i {2^/az) Aa{z, b) + (2\/az) Ba{z, b) 


and 


m 


= r(l+^a-5)e^- ^ 


r(a) 


where 


Ib-i (2\/az) Aaiz, b) - (2v/az) Ba(z, 6)^ , 


A ( u\ 0‘k(z) I bk(z) 

A(2,5)-^-^, Ba(z,b) a-)>oo, 

k—0 k—[) 


(3.1) 


(3.2) 


(3.3) 


7 












valid for bounded values of \z\ and |5|. The coefficients of these expansions can 
be computed by a simple scheme. 

Let 



9{s) 



(3.4) 


The function / is analytic in the strip |3s| < 27r and it can be expanded for 
|s| < 27r into 

CO 

f{z,s) = '^Ck{z)s'^. (3.5) 

k^O 

Then, the coefficients ak{z) and bk{z) of (|3.3|1 are given b£l 


ak{z) = 


m—0 

k 


bk{z) = ^ 


m—0 


{jYl T 2 b)k—mZ Ck-\-m+l{^z) . 


(3.6) 


3.1 Slater’s results for large a 

Slater’s asymptotic expansions of the Kummer functions for large a are given 
in [H §4.6.1], and are in terms of the large parameter a written in the form 


where m > 0 if a and h are real with a > ^b. For programming purposes we have 
simple recurrence relations for the coefficients in the asymptotic expansions. 
Slater’s results for large a are written in the form 




a 2 
6 ’ ^ 




T ! Bk{z) 

zh-i{uz) 2^ + -h{uz) 




an< 


41 


e-i^"z^U {a,b, z^) 


b..l-b 


,2k 




2^U 


r(l a — b) 
Ak[z) z 

9k 

U 


zKb-i{uz) ^ - ^Kbiuz) ^ 


Bk{z) 






The coefficients are given by Aq = 1 and 


,2k 


dt. 


Bk{z) = ^ 

Ak+i{z) = ^^Bk - \B'k{z) + i y eBk(t) dt + Kk, 


(3.8) 


(3.9) 


(3.10) 


^Corrected form of | 15l Eq. (10.3.41)]. 

^Slater’s result for the t/-function contains an error; see |14| and |15l §10.3.6]. 












and Kk is chosen so that Ak+i{z) —?> 0 as z —>• 0. In fact, 


Ao{z) = 1, Bo{z) = 

A,{z) = lib-2)z^ + ±z^ 

Bi{z) = -|6(6 - 2)z - ^z^ + j^z^, 

Mz) = -^{5b - 12){b + 2)z^ + ^{5b - 52)z8 + 
B 2 {z) = ^{5b-12){b + 2){b+l)z^ - 

3^(1756^ - 3506 - 1896)z^ - ^z^^ + 


4 Convergent expansions based on the asymp¬ 
totic representation 

We consider again the representation given in (j3.1l) , but instead of using asymp¬ 
totic expansions of Aa(z, b) and Ba{z, b) in negative powers of a, we expand 


Aa{z, 6) = ^ ak{a, b)z^, Ba{z, 6) = ^ Pk{a, b)z^, 
k=0 k=0 

where the coefficients follow from the recurrence equations (fc > 1) 
f Q!fc-i = 26afe -|- 4(A: -I- l)(fc + b)ak+i — 4(2fc -|- l)/3fc) 

\ /3fc-i = ‘2'bj3k + 4(A: -I- l)(fc -\-2 — b)Pk+i ~ 8a{k + l)afc+i. 
These follow from the system of differential equations 

r 4zy4" -h 46^' - 8zB' + (26 - z)A - 4B = 0, 

\ 4zB" -f 4(2 - b)B' - 8aA' + (26 - z)B = 0, 


(4.1) 


(4.2) 


(4.3) 


which can be obtained by applying Rummer’s differential equation to the rep¬ 
resentations in (EH) and (EJ). 

By solving EH and (13.21) for the A- and i?-functions, and replacing the 
R-Bessel functions by /-Bessel functions, we find that Aa{z, 6) and Ba{z, 6) are 
entire functions of z. Hence, the series in converge for all finite z. 

First values are 


ao 


Po 


a^-^T{a) 

r(a-/l-6)’ 

a{ao - 1) 


0^0 (6^ — 6 -I- 2a) — 2a 
26(1 - 6) ’ 
a (ao (4a — 26 -|- 6^) — 4a -I- 6^) 
26(6-l)(6-2) 


(4.4) 


These values are well defined for 6 = 0,1, 2. 
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For large a we can use [13 Eq. (6.5.72)] for uq and we see that /3o = 1 + 
0{l/a). From a few Maple experiments, we see that a set of other coefficients 
in the expansions in (ED are 0(1/a), and it is likely that we should not use 
forward recursion, when a is large. 

We can change the mixed recursions in (14.21) into two recursion for ak and 
Pk- We have (k > 3) 


p-^ak-s +p- 20 'k -2 +p-iQ!fc-i +poak +Piotk+i = 0, 
q^sPks + q-2l3k-2 + q-iPk-l + QoPk + qiPk+l = 0, 


(4.5) 


where 

P-3 = -(2fc - 1)(2A: + 1), 

P-2 = 8b{2k + l){k — 1 ), 

p_i = —8 + 46 + 8k^ + 24:k — 64A:^ + 32A:‘^ + 126^ — 16fc6 + 166fc^ — 16A:^6^ + 16fc6^, 
Po = 16k{2k — 3)(8a/c^ — 2bk^ — 6^ — 2a + 6), 

Pi = 16k{2k — l)(2k — 3)(fc + 1)(A: + b){—k —1 + 6), 

(4.6) 


q -2 = -2b{2k-l), 

g_i =4(fc-l)(-2fc2 + 62), (4.7) 

go = —8k{2k + 1)(A: — l)(4a — 6), 

gi = 16fc(fc — l)(fc + l)(fc + b){k + 2 — 6). 

In order to decide which of the directions of recursion is stable (if any), it is 
important to study the possible existence of minimal solutions. A simple tool 
that can be used is the Perron-Kreuser theorem [3 (see for instance [5] for a 
more recent account of this important result), which in some cases gives simple 
answers to this question. This is not the case here, however. 

To see this, first we notice that the coefficients of the recurrence behave as 
k —>■ +00 as follows: 

gi ^ 16fc^, go ~ —16(4a — 6)fc^, g_i ^ —8fc^, g _2 —46fc, g _3 ^ k. (4.8) 

The coefficients pi have similar behaviour, and we have pi ^ —‘ikqi as k ^ +c». 
Considering the Newton-Pusieux diagram corresponding to this asymptotic be¬ 
haviour, we conclude that all the solutions Pk of the recurrence relation satisfy: 

limsup(/c!|/3fc|)i/'= = i, (4.9) 

k—¥+oo ^ 

where the constant \ is the absolute value of the (degenerate) solutions of the 
characteristic equation 16A^ — 8A2 + 1 = 0, and the coefficients in this equation 
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Figure 1: Relative errors in the computation of U{a, b, z) using the expansions 
of (14.11) with a fixed number (TV = 20) of ak and f3k coefficients for computing 
Aa{z,b), Ba{z,b) in (|3^ . 


come from the coefficients gi, q-i and q-^. The same result is true for ak- This 
is all the information the Perron-Kreuser theorem gives in this case; it is not 
possible to infer the existence of minimal solutions only with this result. 

A scaling of the sequences ak and /3fc in order to minimize the risk of over¬ 
flow/underflow in the recurrence can be used, by defining ak = k\2^ak and 
similarly for the /3fc’s. Numerical experiments show that for small a, the for¬ 
ward computation appears to be stable while for large a this is no longer true. 

Although, the Perron-Kreuser theorem does not give information on the ex¬ 
istence of minimal solutions in this case, it is possible to test numerically the 
existence of minimal solution. One can start the recurrence with large values of 
k and arbitrarily chosen starting values, apply the recurrence in the backward 
direction and check whether the same value ai/ao (or l3i/f3o) is obtained irre¬ 
spectively of the initial values. If this is the case, this means that a minimal 
solution must exist. We have checked this fact for several sets of starting values 
and values of the parameters and it seems that such a minimal solution exists; 
however, it does not appear to be the requested solution (because the values do 
not correspond to (14.41) 1. 

Therefore, backward recurrence does not seem to be of help. Contrarily, 
there is a range of variables for which the the forward recurrence yields correct 
results for the expansions. For instance, when \z\ < 1, |6| < 1 (where the b- 
values should be neither too close to 0 nor to 1) we obtain an accuracy close to 
double precision (or better) for values of the product \za\ smaller than 10 for 
computing the U function, as illustrated in Figure [TJ In this figure we show, 
for three values of the b parameter, the relative errors in the computation of 
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Figure 2: Number of terms needed in the expansions of 63) to compute 
U{a,b,z) with an accuracy better than The value of the parameter b 

has been hxed to 0.4. 


U (a, b, z) using the expansions of (14.111 with a fixed number of terms (iV = 20) for 
computing Aa{z,b), Ba{z,b) in (14.ip . On the other hand, when the demanded 
accuracy in the computation of U{a,b,z) is fixed to 10“^“^, for example, the 
number of terms needed in 63 seems not to be larger than 10, as shown 
in Figure [2] for & = 0.4; this also applies for other values of the b parameter, 
excluding values very close to 0 or 1. 
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